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In Einstein-aether theory and Horava gravity, a timelike unit vector is coupled to the spacetime 
metric. It has previously been shown that in an exponentially expanding homogeneous, isotropic 
background, small perturbations of the vector relax back to the isotropic frame. Here we investigate 
large deviations from isotropy, maintaining homogeneity. We find that, for generic values of the 
coupling constants, the aether and metric relax to the isotropic configuration if the initial aether 
hyperbolic boost angle and its time derivative in units of the cosmological constant are less than 
something of order unity. For larger angles or angle derivatives, the behavior is strongly dependent 
on the values of the coupling constants. Generally there is runaway behavior, in which the anisotropy 
increases with time, and/or singularities occur. 
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I. INTRODUCTION 

When the phenomenology of theories with a preferred 
frame is studied, it is generally assumed that this frame 
coincides, at least roughly, with the cosmological rest 
frame defined by the Hubble expansion of the universe. 
Observations place strong bounds on frame dependent 
effects which would presumably grow with the relative 
velocity of the preferred frame and the velocity of the 
earth (which moves at the "low" speed of 10~ 3 c rela- 
tive to the Hubble frame). In a particular theory with 
preferred frame effects, the dynamical alignment of the 
frame (or frames) can be studied to determine stability 
of cosmic alignment, as well as to characterize the range 
of initial conditions that could be expected to naturally 
align. 

In this paper we examine this question in the case 
of Einstein-aether theory[T] and in the IR limit of (ex- 
tended) Horava gravity [21 [3]. Einstein-aether theory just 
consists of general relativity coupled, at second derivative 
order, to a dynamical timelike unit vector field u a , the 
aether. In Horava gravity, the aether vector is assumed 
to be hypersurface-orthogonal, i.e. it is the unit normal 
to level sets of a scalar time function. Various forms of 
"Horava gravity" have been discussed in the literature. 
Here we refer exclusively to the one related to Einstein- 
aether theory as just explained. (This corresponds to 
the so-called "non-projectable" version, where the lapse 
function N is an arbitrary function of spacetime, and 
includes in the Lagrangian a term proportional to the 
square of the gradient of IniV.) Every hypersurface- 
orthogonal Einstein-aether solution is a Horava solution. 
All the solutions to be considered in this paper are of this 
type. 

The alignment of the aether has been studied before 
in the context of linearized perturbations. The question 
was first addressed, indirectly, by Lim[3], who found that 
all perturbations of the aether decay exponentially in a 
de Sitter background. In particular, this result applies 
to the homogeneous modes. Subsequent workj5] [6] con- 
firmed this result, but in Ref. [B] it was found that under 



some circumstances, after inflation, velocity perturba- 
tions might grow to be "mildly relativistic" and could still 
possibly be compatible with observations. In all these 
analyses, it is assumed that the aether is aligned in a 
background solution, and the behavior of perturbations 
is studied. 

Kanno and Soda (KS) approached the question from a 
different point of view. In the Appendix of Ref. [7] they 
examined homogeneous but anisotropic solutions in the 
presence of a positive cosmological constant, with three 
orthogonal principal directions of expansion, and with 
the aether tilted in one of the principal directions. [This 
corresponds to Bianchi type I (Kasner-like) symmetry] 
They showed that, to linear order in the anisotropy, the 
system relaxes exponentially to the isotropic, de Sitter 
solution. Since their analysis was carried out just to lin- 
ear order in the anisotropy, it is in fact just a special case 
of the above-mentioned perturbative treatments. 

In this paper we adopt precisely the setting of the KS 
analysis but we include the full nonlinear dynamics. We 
characterize the range of initial data that relax to an 
isotropic solution. Generically, the aether aligns pro- 
vided the initial boost angle and its time derivative in 
units of the cosmological constant are less than something 
of order unity. The precise stability bounds depend on 
the values of the coupling parameters in the Lagrangian 
defining the theory. 



II. BIANCHI TYPE I EINSTEIN-AETHER 
COSMOLOGY 

Einstein-aether theory is general relativity (GR) cou- 
pled to a dynamical timelike unit vector field. In terms of 

the metric g a b of signature (H ) and the unit vector 

field u a it is defined by the action 

S = / d 4 x ^ (i? + 2A + K ab mn V a u m V b u n ) 

(1) 
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where R is the Ricci scalar, A is a cosmological constant, 
and the tensor K ab rnn is given by 
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and ci , . . . , C4 are dimensionless coupling parameters that 
define the theory. Since u m is constrained to be a unit 
vector, the action need only be stationary under vari- 
ations orthogonal to the aether, u m Su m = 0. In this 
paper for simplicity we omit any matter couplings, since 
the cosmological constant suffices to source the overall ex- 
panding solution and it models the conditions that would 
have pertained in an inflationary early universe. It would 
be straightforward to add radiation or matter or some 
form of quintessence to the model. 



A. Bianchi type I symmetry 

Following KS, we specialize to Bianchi type I space- 
times, i.e. to metrics that are homogeneous and spatially 
flat, with three commuting translation symmetries, 



ds 2 = N 2 (t)dt 2 - e 



2a(t) 



-4a 



-(*) 



dy 2 



-Udz% (3) 



We also assume that the aether vector is tilted only in 
the x-direction, 



f 



N(t) 



cosh6>(i) d t + e- a (*) +2<T +(*) sinh 9{t) d x . (4) 



The hyperbolic angle 9 measures the boost of the aether 
relative to the rest frame of the homogeneous, flat spa- 
tial sections, i.e. the "homogeneous frame". The metric 
is determined by four functions. The lapse N(t) specifies 
the flow rate of proper time with respect to t in the ho- 
mogeneous frame. Co-moving lengths L in the x, y, and 
z directions all have different expansion rates, L/L. The 
sum of these is 3d, which is also the fractional rate of 
change of co-moving volume, V/V. The quantity 2v / 3o'_ 
is the difference between the expansion rates in the two 
transverse directions y and z, while 3d+ is the difference 
between the average of these and the rate in the x direc- 
tion. 

The vector field Q is in effect just two dimensional; 
hence, like all such vector fields, it is hypersurface or- 
thogonal. According to the analysis of Ref. [8j, this 
means that the solutions to the field equations discussed 
here are also solutions to the field equations of Horava 
gravity. Hence our results apply to the cosmology of 
that theory as well. The hypersurface orthogonality also 
means[5] that the action is unchanged under c\ — > Oy +6, 
c 3 C3 — t5, and C4 — > C4 — 5, so the system depends 
on these three coupling parameters only through the two 
invariant combinations c± +C4 and c\ +C3. For notational 
compactness we shall make use of these quantities, and 
also drop the subscript 2 on C2- 



c\ + Ci, 



ci + c 3 , 



C2- 



(5) 



[These parameters correspond respectively to the param- 
eters a, j3, A' of the action for Horava gravity, Eq. (5.72) 
in Ref. [TU].] 

When the fields have this symmetry structure, the ac- 
tion takes the form (up to a total derivative) 



S : 
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where q l f-> (9, a, er+, er_). Here and below the time de- 
pendence of the dynamical variables is implicit, the dot 
denotes derivative with respect to i, and indices i, j, . . . 
label the four dynamical variables. The nonzero compo- 
nents of the symmetric array ffy are given by 



Hee 



2(6 + c + (a - b - c) cosh 2 9) 
2 (a — b — 3c) cosh 9 sinh 9 
4(-a + b) cosh 9 sinh 9 
2(-6 - a + (a - 3b - 9c) cosh 2 1 
-4a sinh 2 9 

4(3 - 2a + (2a - 36) cosh 2 9) 
12(1 - 6 cosh 2 9) 



(7) 



The 9 dependence of Hij(9) will also be suppressed. 
The dynamics is symmetric under the inversion {9,9) — > 
(-9,-9). 

A key to the general behavior of solutions is the invcrt- 
ibility of Hij, whose determinant can be written in the 
form 

detH = -1728a(I-6) 2 (2 + 6 + 3c) 
x (vl + (I - vl) cosh 2 9) (vl + (I — v 2 2 ) cosh 2 9) (8) 



with 



1-6' 



(b + c)(2-a) 
a(l-6)(2 + 6 + 3c)' 



(9) 



The constants vq and V2 are, in fact, the speeds of the 
spin-0 and spin-2 modes linearized around flat space (de- 
fined with respect to the background aether frame), in 
both Einstein-aether theory[T] and Horava gravity. Note 
that the vq factor is proportional to H 

It may initially be surprising that the spin-0 and spin-2 
perturbations play any role in a homogeneous cosmology. 
But while the metric (|3| is homogeneous on constant t 
surfaces, it has x dependence on surfaces orthogonal to 
the aether. In fact, the form of ^ can be understood 
from simple kinematic considerations as follows. 

The linearized, diagonalized, action for a mode tp with 
speed v is proportional to gfydaijjdbip , where gfK oc 

u a u b + v 2 (g ab — u a u b ) is the effective metric for that 
mode. Since we consider only homogeneous fields, which 
depend on t alone, the only component that enters is 
gf v) a (v 2 g tt + (l-v 2 )u t u t ) cx (v 2 + (l-v 2 )cosh 2 (9). The 
symmetry and hypersurface orthogonality of the aether 
permit only one spin-2 mode and the spin-0 mode. The 
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spin-2 mode is governed by the shear cr_ , which describes 
the gravitational wave mode transverse to the tilt and 
with polarization aligned with the y and z symmetry 
axes. This explains the form of the determinant 

If the coupling constants (a, b, c) are such that one of 
the mode speeds exceeds unity, then there is a value of 
9 for which the determinant of Hij vanishes, correspond- 
ing to the condition g^-, = 0. When the aether reaches 
this hyperbolic tilt angle, the propagation cone of that 
mode becomes tangent to the constant t surface, so that 
becomes a valid constant phase surface for the mode. 
In other words, the mode propagates instantaneously on 
the constant t surface. Beyond this aether tilt the kinetic 
energy of the mode becomes negative, so the system is 
unstable. In section |TV| we shall discuss the implications 
of this phenomenon for the cosmological dynamics. 



B. Equations of motion for A = 

Although our main interest is in the case of exponential 
expansion driven by a positive cosmological constant, we 
begin by looking first at the simpler case of vanishing 
A. Then variation with respect to the lapse N yields the 
initial value constraint 



0. 



(10) 



As usual in general relativity, if this constraint is satisfied 
at one time, then the rest of the equations of motion 
imply that it remains satisfied for all time. We can choose 
the nontrivial lapse 



N 



,,3o 



(11) 



to eliminate the a dependence in the action. Then 
the dynamics becomes that of affinely parametrized null 
geodesies on the configuration space (6, a, a +1 a J) with 
respect to the metric that depends only on 9. 
It is convenient to define momenta by 



which can be solved for the velocities, 



(12) 



q l =H 
when the inverse H %0 of H 



menta, the constraint ( 10 1 reads 



Pi, (13) 
exists. In terms of the mo- 



H l3 ViPi = 0. 



(14) 



Since depends only on 9, the momenta p a and p± are 



conserved. Moreover, the constraint (14) is a quadratic 
equation in pg that can be solved for pe(9;p a ,p±) (there 
are generically two roots or none). Having solved three 
of the four evolution equations, as well as the constraint 
equation, the fourth evolution equation, for pg, is now 
redundant. The dynamics for 9 is thus reduced to a first 
order differential equation, 



9 = H' 



k 



Pk 



F(9; Pa ,p ± ). 



(15) 



Once the evolution of 9 is known, the remaining variables 
a and cr± are determined by integration of the first order 
equation (Il3|) . The character of the evolution of 9 can be 



seen by inspection of a plot of the graph of the function 
F defined in ((TBI). 



C. Equations of motion for A ^ 

For nonvanishing A, the variation with respect to the 
lapse N yields the initial value constraint, 



//,,</''/' + 4A = 0. 



(16) 



Because of the A term in the action, the lapse ( 11 ) is no 



longer the most convenient, and it is simpler to just use 
N = l. (17) 
The Euler-Lagrange equations with this gauge choice are 

d (e 3a H ijq i) ~ diie^^HuqKq 1 - 2A)) - 0. (18) 



dt 



The individual components i = 6, a, ± read 



e"~H aj q 3 ) 



dt 
d 
dt 



V"' //, ,.,,,/'</ 
-12 V 

0, 



i l n3 



(19) 
(20) 
(21) 



where in ( 20 ) the constraint ( 16 1 was used 



Again, it is sometimes convenient to express the field 
equations in terms of the "momenta" (12 1. (These are 
not precisely the conjugate momenta anymore since the 
factor e 3a is not included, but we will nevertheless refer 



to them as momenta.) Then the constraint (16) takes 
the form 



H v p iPj + 4 A = 0, 

and the equations of motion become 

pe = -3ape - \H 13 g piPj 
p a = —'3ap a — 12A 
p ± = -3ap ± 



(22) 

(23) 
(24) 
(25) 



Using ( 13 ) we can express a in terms of 9 and the mo- 
menta, 



a = H ak Pk . 



(26) 



Moreover, we can solve the quadratic constraint equa- 
tion ( |22[ ) for p a , so that a and p a can be eliminated 
completely from the dynamical system. 



III. LIMITING CASES 

In this section we discuss various special cases and lim- 
its of the theory. 
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A. General relativity 



the Hij array reduces to the nonzero elements 



If we reduce to the pure GR case, 9 is not present, and 
a = b = c = 0. Then Hj is diagonal and constant, with 

H ++ = H = -H aa = 12. The constraint ((22]) then 

becomes 



2,2,2 
P a +P++P- 



-48A. 



(27) 



The only isotropic (p± = 0) solutions are Minkowski 
spacetime (with A = 0) and de Sitter spacetime (with 
a = J A/3, A > 0). In the anisotropic case, A = yields 
the Kasner solutions, and A ^ yields a generalization 
of those. 



B. 9 = solutions 

Next we characterize the solutions in which 9 = for 
all times, i.e. in which the aether remains orthogonal to 
the constant t homogeneity surfaces. Then, although the 
aether has no dynamics, its couplings in the action ([!]) 
still contribute to the field equations and we therefore 
have something different from GR. 

There are no terms linear in 9 or 9 alone in the ac- 
tion (J6] , and terms of quadratic or higher order in these 
quantities will obviously not contribute to the equations 
of motion if 9 = 9 = 0. If this condition holds initially it 
is therefore preserved for all time, and for characterizing 
these solutions it is consistent to simply set 9 = in the 
action. Then H^ is diagonal and constant, the relevant 
components being 



H aa = -6(2 + 6 + 3c) 
H ++ = 12(1-6) 
ii+_ = 12(1-6). 



(28) 



In the isotropic case o~± = 0, the system is then equiv- 
alent to GR with a rescaled cosmological value of New- 
ton's constant [TQ H2], G cosmo = G/(l + (6 + 3c)/ 2), 
and with A replaced by A' = A/(l + (6 + 3c)/2). This 
isotropic solution is the spatially flat slicing of de Sit- 
ter spacetime with Hubble constant H = \/A'/3. The 
aether becomes singular because of infinite stretching on 
the past horizon. In the anisotropic case, interestingly, 
there is no equally simple relation to GR: the presence 
of the isotropic aether induces different rescalings of the 
kinetic energy associated with expansion and shear. 



C. Linearized anisotropy 

If one drops all terms in the action ^ of higher than 
quadratic order in the anisotropic coordinates 9 and <r±, 



Haa 

H++ = 12(1-6) 
H— = 12(1 - 6). 



2a 

2(a-b-3c)d 

-6(2 + 6 + 3c) + 2{a - 36 - 9c) 



(29) 



Keeping only linear order terms in the anisotropy, the 
equation of motion for 9 reduces to 



9 + 3a9 + 2a 2 9 = 0. 



(30) 



To zeroth order in the anisotropy the solution to the con- 



straint ( 16 1 is 



A 



3[l + (6 + 3c)/2] 



(31) 



and a in (30) can be replaced by this value. Then (30) 



is the equation found by KS [7J. They pointed out that 
the coefficient of 9 is positive provided the effective grav- 
itational coupling is positive, in which case this is the 
equation of a damped harmonic oscillator. In fact, the 
oscillator is overdamped, with eigenmode decay rates a 
and 2a. This implies that 9 relaxes to zero as the universe 
expands. 



D. p± = solutions 



The equation of motion ( 25 ) is solved for p± by 



P±(t) = e~ 3Q(t) P± ,o, 



(32) 



where p± t o is an integration constant. Using this, the 
remaining equations of motion ( 23|24 1 involve only the 
variables {9,pg, a, a,p a ). Moreover, as mentioned at the 



end of section II C Eq. (26) can be used to eliminate 
a, and the constraint can be solved for p a . However 
this requires that the function a(t) be determined by the 
previous values of (9,pe,a) via J dif a(t'), which does 
not yield an evolution equation that is local in time. 

If the co-moving volume is expanding then, according 
to (32), p±(t) is driven to zero. It is therefore a useful 
limiting case to set p± = from the beginning. Then 
there is no remaining a(t) dependence, and the system 
can be reduced to the 9 degree of freedom alone. We now 
explain in detail how this is achieved. 

We assume now that p± = 0. Since H~ k = except 
for k = — , it follows from (13) that b~ = 0. Thus, the 



two transverse dimensions must have the same expansion 
rates. In contrast, it does not follow that <r + = 0, since 
we have in this case 



H 



Pe 



H^ 



Pa 



(33) 



As might be expected, 
ric to be isotropic (c 



it is inconsistent for the met- 
= 0) when the aether is tilted 
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(0^0). (However, note that & + is of second order in the 
9 anisotropy.) The expansion rate in the tilt direction 
must generally differ from that in the transverse direc- 
tion. Eliminating <j + via (33), the system reduces to the 
variables (9, pg, a, p a ). 

One can write this system in terms of the velocities, 
i.e. in t erm s of the variables (9, 9, a, a), by using the con- 
straint ( 16 1 to solve for a + in terms of (9, 9, a), and sub- 



stituting that into the Euler-Lagrange equation ( 18 ). But 
for the purpose of making a (9, 9) phase portrait of the 
evolution, it appears more neat to organize the equations 
as follows. 

The idea is to solve for the momenta (pg,p a ) in terms 
of the velocities (9, a), and then to use the equations 
that were expressed in terms of momenta. To this end, 
we introduce capital indices A, £?, ... to refer to the two 
coordinates 9 and a, we define the contravariant tensor 
h AB to be the restriction of , 



AB 



H 



AB 



(34) 



and we denote by Hab the inverse of h AB . Then from 
(131 we have 



q A = h^p Bl 
which can be inverted to yield 

Pa = h A Bq B ■ 



Using ( 36 1 the constraint ( 22 1 becomes 



h AB q A q a + 4A = 0, 

which can be solved as a quadratic equation for d( 
thus eliminating a. Explicitly, we have 



(35) 



(36) 



(37) 



-hg a 9 ± J(hg a 9) 2 - h aa {hgg9 2 + 4A) 

a= ¥ . (38) 

There are generically two solutions or no solutions. If 
h aa {hggt r+4A) < there are two solutions, one in which 
the volume is expanding (d > 0) and the other in which 
it is contracting (a < 0). Note that for 9 = the solu- 
tions are simply a = ±y/—4A/h aa . As h aa — > CP this 
diverges, and no solution exists when h aa > 0. 

The remaining task is to find an equation for 9. For 
this purpose we can use (13) to write 9 = h eB pB] hence, 



>Pb 



Pb- 



(39) 



Then using (36) and the equations of motion ( 23|24 | we 
find, after some manipulation, 

9 = -3a9 -12Ah ea - h eA h A B,eq B (40) 

(41) 



n* 1 n-AB.eq q 



Together with ( 38 ) this yields a dynamics reduced to just 



the 9 degree of freedom, which can be visualized in a 
phase portrait. 



IV. GENERIC BEHAVIOR 

We have seen that, in a small enough neighborhood of 
9 = 0, the dynamics relaxes to the 9 = case, provided 
the values of a, 6, and c are such that the effective gravi- 
tational coupling constant is positive. On the other hand, 
once 9 is sufficiently large, the character of the dynamical 
system can obviously change dramatically because of the 
growth of the hyperbolic trigonometric functions in the 
components of Hij Q . 

A general feature mentioned earlier is that det H ^ 
vanishes if either the spin-0 or spin-2 propagation cone is 
tangent to the constant t surface. The conditions deter- 
mining these angles can be expressed as 



coth#o — Vq, COtll#2 = V2 



(42) 



where t>o,2 are the mode speeds (9b. At either of these 
ang les is not invertible, so the equation of motion 
(18) cannot be solved for if. As such a value of 9 is 
approached, at least one second derivative component 
would generally diverge. Hence generically there can be 
no smooth evolution across the degenerate values of 9. 
The dynamics may run into a singularity there, or it may 
"bounce" before reaching such a value of 9. 

There is a solution 9* to each of the equations in ( [42] ) 
as long as the corresponding mode speed v is greater than 
unity. The larger of the mode speeds defines the smaller 
of the critical angles. The critical angles are of order 
unity unless v is either very large or very close to 1. In 
these limits we have 



l/v 
-|ln(u- 1) 



for 
for 



v » 1 
v - K 1 



(43) 



In particular, the degenerate value 9i is real only if < 
b < 1, and is of order unity unless b is very close to either 
or 1. For instance, for b = 0.01, 0.9, or 0.99, we have 
6*2 — 3, 0.3, or 0.1 respectively. The degenerate value 9q 
is of order unity for generic values of a, b, c with no large 
hierarchy amongst them. [If a, 6, c are all much smaller 
than 1, then vq ~ (b + c)/a, so 9q sa coth _1 ((6 + c)/a).] 
We infer that exotic behavior, including singularities or 
runaway solutions, may typically occur for aether boost 
angles of order unity. 



A. Restriction to physically viable couplings 

There are three independent coupling constants that 
affect the solutions we are studying in cither Einstein- 
aether theory or Horava gravity, but stability and obser- 
vational constraints restrict the range of physically viable 
values. 



1. Einstein- aether couplings 

As summarized in Ref. [T], and C4 should be de- 
termined by ci and C3 such that the preferred frame 
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parametrized post-Newtonian (PPN) parameters ai.2 
vanish (or are very small compared to unity). More- 
over, when ai.2 vanish, stability and positive energy of 
linearized modes and the absence of vacuum Cherenkov 
radiation by ultra high energy cosmic rays require < 
c + < 1 and < c_ < c+/(3(l — c + )), where 



c± = ci ± c 3 . 



(44) 



(In terms of a, 6, c these conditions correspond to < 
6 < 1, < a < 26/(4 - 36), and c = -(a + 6)/3, or, 
equivalents, < 6 < 1, 6(6 - 2)/(4 - 36) < c < -6/3, 
and a = — b — 3c.) We shall label the examples by their 
c± values. In particular, the vacuum Cherenkov con- 
straint requires that all the mode speeds be greater than 
or equal to unity, so except in the case where they are 
exactly unity, there are values of the tilt angle where the 
dynamics is degenerate. The degeneracy at 02 is relevant 
to the dynamics only if <j_ =/= 0. 

The remaining observational constraint one might ap- 
ply is that the radiation damping rate for a binary pul- 
sar system agrees with the rate in GR, which agrees with 
observations within the present relative uncertainty of 
about 0.002. The results of Ref. [13] establish that this 
constraint is satisfied for generic small values Cj < 0.001, 
and if Cj < 0.01 — 0.1 it is satisfied if c_ « 0.18 c + . For 
these values the spin-0 mode speed is v ~ 1-36 and the 
critical boost angle is 9q ~ 0.94. 



2. Horava gravity couplings 

The constraints on the couplings in Horava gravity are 
the same as in Einstein-aether theory except for the PPN 
constraints = 0, which now are equivalent to a — 
26 [10] . The other constraints are < 6 < 1 and (6 + 
c)/(6(2 + 6 + 3c)) > 1. Given the first of these, the second 
is satisfied in two regions: (I) c > (6 + 6 2 )/(l — 36), 
6 < 1/3, and (II) c < -2/3 - 6/3, c> (6 + 6 2 )/(l - 36) 
when 6 > 1/3. The radiation damping rate has not yet 
been calculated in the Horava case. 



B. Phase portraits 



As discussed in Sec. Ill D| in the case when p± = 
one can reduce the dynamics to the 9 degree of freedom. 
Then the dynamics can be displayed as a phase portrait 
in the (9, 9) plane, exhibiting the flow of the vector field 
(9, 9). This serves to illustrate the general features of the 
dynamics discussed above. 

Two examples are shown in Figs. [TJand[2] In both of 
these the parameters C2 and C4 are chosen so as to sat- 
isfy the PPN constraints of Einstein-aether theory, and 
c± satisfy the remaining constraints other than that of 
gravitational radiation damping. Figure[TJis qualitatively 
similar to the phase portrait for the case c+ = 1/10, 
c_ = 0.18c + , which is at least close to satisfying the ra- 
diation damping constraint. It is also similar to the case 
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Figure 1: Stream plot of the vector field (6,6) on the (6,6) 
plane, with c + = 1/10 and c_ = 1/40, and p± = 0. In this 

± 1.16. 




Figure 2: Stream plot with c+ = 1/2 and c_ = 1/4, and 
p± = 0. In this case do ~ 1.3, but already for 6 > 1.1 and 
8 — there is no solution to the constraint equation for a. 



(a,b,c) = (2/10,1/10,3/10) which satisfies the Horava 
constraints other than the (unknown) radiation damping 
one. 

In the case illustrated by Fig. [TJ runaway behavior 
occurs if 9 or 9 is sufficiently large. Numerical evo- 
lution of this example suggests that some of the flow 
lines end at curvature singularities. The simplest space- 
time scalar in this setting is the expansion of the aether, 
V ' a u a = 3d cosh 9 + ^sinh^. For example, with initial 
data 9 = and 9 = 4, the evolution runs away to large 
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9. For another example, with initial data 9 — 1.5 and 
9 ~ 1.8, the evolution runs away to large a. 

The singular behavior seen in these solutions might be 
related to what is seen in some homogeneous anisotropic 
cosmologies with tilted perfect fluid matter (and vanish- 
ing cosmological constant) [14] . It may appear inconsis- 
tent with the cosmological "no-hair" theorem proved by 
WaldfTB]. which showed that in the presence of a positive 
cosmological constant A, all expanding Bianchi-type cos- 
mologies (except Type IX) evolve toward the de Sitter 
solution with time scale y/3/A. But that result assumed 
that the dominant and strong energy conditions hold for 
the matter stress tensor. These conditions do not gener- 
ally hold for the stress tensor associated with the aether 
part of the action ([!]). 

V. CONCLUSION 

The question driving this investigation was whether it 
is natural for the aether to be aligned with the isotropic 
frame of a homogeneous, isotropic cosmology in Einstein- 
aether theory or Horava gravity? We addressed this ques- 
tion by studying the dynamics of a tilted aether in a 
homogeneous anisotropic Bianchi type I cosmology with 
a cosmological constant. We found that generically the 
aether does align provided its tilt angle and the time 
derivative of its tilt angle in units of the cosmological 



constant are smaller than something of order unity. This 
extends the linearized stability result of KS 17] to a fi- 
nite basin of attraction whose precise shape depends on 
the coupling parameters of the theory, and in some cases 
the basin appears to be much broader than order unity. 
Outside of this basin, the solutions exhibit runaway or 
singular behavior of one or more of the variables. Some 
of this behavior occurs when the propagation cone of ei- 
ther the spin-0 or spin-2 mode is tilted enough to meet 
the homogeneous constant t surface. We do not know 
whether similar behavior would persist if the homoge- 
neous symmetry condition were dropped. 

Our findings show that the fate of a universe with 
Bianchi I symmetry depends heavily on the initial tilt of 
the aether. Perhaps the question of the initial tilt could 
be addressed from the standpoint of quantum cosmology, 
for example via the "wave function of the universe" or via 
the distribution of initial conditions for chaotic inflation. 
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